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Abstract 

Molecular dynamics (MD) simulations have been carried out to investigate 
the slip of fluid in the lid driven cavity flow where the no-slip boundary con- 
dition causes unphysical stress divergence. The MD results not only show the 
existence of fluid slip but also verify the validity of the Navier slip boundary 
condition. To better understand the fluid slip in this problem, a continuum 
hydrodynamic model has been formulated based upon the MD verification of 
the Navier boundary condition and the Newtonian stress. Our model has no 
adjustable parameter because all the material parameters (density, viscosity, 
and slip length) are directly determined from MD simulations. Steady-state 
velocity fields from continuum calculations are in quantitative agreement with 
those from MD simulations, from the molecular-scale structure to the global 
flow. The main discovery is as follows. In the immediate vicinity of the cor- 
ners where moving and fixed solid surfaces intersect, there is a core partial-slip 
region where the slippage is large at the moving solid surface and decays away 
from the intersection quickly. In particular, the structure of this core region is 
nearly independent of the system size. On the other hand, for sufficiently large 
system, an additional partial-slip region appears where the slippage varies as 
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1/r with r denoting the distance from the corner along the moving solid sur- 
face. The existence of this wide power-law region is in accordance with the 

asymptotic 1 jr variation of stress and the Navier boundary condition. 
47.11. +j, 68.08.-p, 83.10.Mj, 83.10.Ff 
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I. INTRODUCTION 



A crucial ingredient in the continuum hydrodynamics is the boundary condition of fluid 
flow past a solid surface. The no-slip boundary condition, i.e., zero relative velocity between 
the fluid and solid at the interface, is a core concept in fluid mechanics [1]. In molecular 
dynamics (MD) simulations, however, a small amount of relative slip between the fluid 
and the solid surface is generally detected [2-5]. Such slip can be accounted for by the 
Navier boundary condition (NBC), whereby the slip velocity is proportional to the tangential 
viscous stress and the degree of slip is measured by a slip length [2-5]. As the relative slip 
is extremely small in macroscopic flows, the NBC is practically indistinguishable from the 
no-slip boundary condition in most situations. 

When applied to the immiscible two-phase flow of moving contact line, where a fluid-fluid 
interface intersects the solid wall, the no-slip boundary condition would cause non-integrable 
diverging stress and unphysical infinite dissipation, which directly imply the breakdown of 
the no-slip boundary condition [6]. In the past two decades, MD simulations have shown 
fluid slip in the molecular-scale vicinity of the moving contact line [7,8]. Recent evidences 
have shown the slip velocity profile obtained from MD simulations to be accountable by the 
generalized Navier boundary condition [9], in which the slip velocity is proportional to the 
total tangential stress — the sum of the viscous stress and the uncompensated Young stress; 
the latter arises from the deviation of the fluid-fluid interface from its static configuration. 

The no-slip boundary condition also runs into trouble when applied to the driven cavity 
flow, where a rigid plane slides steadily over another, with a constant inclination angle 
[1,10]. This geometry readily shows that the no-slip boundary condition would cause the 
same non-integrable diverging stress as in the problem of moving contact line. Near the 
corner where fixed and moving solid surfaces intersect, the velocity variation becomes very 
fast because two different velocities are assumed at the two solid surfaces. Moreover, a 
velocity discontinuity occurs at the corner if the no-slip boundary condition is everywhere 
applied. This is the origin of the non-integrable 1/r stress at r — > 0. A slip boundary 
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condition is therefore imperative. 

Koplik and Banavar first used MD simulations to explore the small-scale structure of 
the driven cavity flow [10]. Their results indicate that slip occurs in the corner region. To 
uncover the slip mechanism, they measured the microscopic tangential stress at the fluid- 
solid interface. This stress measurement led to the conclusions that local non-Newtonian 
region exists in a low-shear and otherwise Newtonian flow (at low Reynolds number) and 
that the NBC is not valid. 

The purpose of this paper is to uncover the slip boundary condition and formulate a 
continuum hydrodynamic model for the driven cavity flow, from which we will answer the 
intriguing question: In a mesoscopic or macroscopic system, what is the slip profile which 
consistently interpolates between the inevitable slippage in the immediate vicinity of the 
corner and the no-slip boundary condition that must hold at mesoscopic/macroscopic regions 
far away? Success would lead to new understanding to slip and dissipation in restricted 
geometries [11]. 

The continuum hydrodynamic modeling requires a slip boundary condition and a mo- 
mentum transport equation. We have carried out MD simulations similar to those performed 
for the moving contact line problem [9]. In contrast to the conclusions in Ref. [10], the NBC 
is found to be governing the slip of fluid relative to the solid. The shear stress is verified to 
be Newtonian in the vicinity of the corner where velocity and stress variations are extremely 
large. Technically, we denote a molecular boundary layer of fluid at the fluid-solid interface. 
We then perform stress measurement using a method that is reliable near the interface [9]. 
Velocity and stress data collected at the boundary layer provide molecular evidence for the 
validity of the NBC. We emphasize that unlike the Couette flow simulated in Ref. [3-5], 
here the driven cavity flow requires a local verification of the NBC, because the slip velocity 
varies along the solid surface and this variation becomes very fast in the corner region. As 
for the Newtonian behavior, local velocity and shear stress are measured everywhere in the 
fluid, showing that the shear stress is proportional to the local shear rate. 

Based upon the conclusions drawn directly from the simulated cavity flow, a continuum 
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hydrodynamic model has been formulated, comprising the Navier-Stokes equation and the 
NBC. With all the material parameters directly determined from MD simulations, our model 
has no adjustable parameter. Numerical calculations have been carried out to produce 
continuum results for comparison with MD results. It is shown that in a wide range of 
Reynolds number, MD and continuum flow fields agree well, from the molecular-scale corner 
region (with large slip) to the large-scale outer region (with vanishing slip in a large system). 
The largest Reynolds number ever reached is « 50, at which deviation from the Stokes flow 
is clearly noted (see Sec. VI D). 

The paper is organized as follows. We first describe the details of the MD simulations in 
Sec. II. We then outline in Sec. Ill our MD approach to the verification of the Navier slip 
boundary condition. The MD results, later used for hydrodynamic modeling, are presented 
in Sec. IV. A continuum hydrodynamic model is formulated in Sec. V. The numerical 
algorithm is also briefly described. In Sec. VI there is a systematic comparison of the MD 
and continuum hydrodynamics results. The paper is concluded in Sec. VII with a few 
remarks. 

II. MOLECULAR DYNAMICS SIMULATIONS 

The purpose of carrying out MD simulations is threefold: (1) To uncover the boundary 
condition governing the driven cavity flow (Sees. Ill and IV); (2) To determine the material 
parameters (e.g., viscosity and slip length) in our hydrodynamic model (Sec. V); (3) To 
produce flow fields for comparison with the continuum hydrodynamic solutions (Sec. VI). 

We consider a single fluid confined in a two-dimensional (2D) cavity formed by two 
horizontal walls in the xy plane and two vertical walls in the yz plane (see Fig. 1) [10]. The 
cavity measures L along x and H along z, with the periodic boundary condition applied 
along y. The fluid is sheared by moving the upper and lower walls with the same speed 
V w along the ±x directions, respectively. Each of the four walls is constructed by two to 
four [001] planes of an fee lattice, with each wall molecule attached to the lattice site by 
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a harmonic spring. The mean-squared displacement of wall molecules is controlled to obey 
the Lindemann criterion. Interaction between the fluid molecules separated by a distance r 
is modeled by a Lennard- Jones (LJ) potential 



where e and a are the energy scale and range of interaction, respectively. The wall-fluid 
interaction is modeled by a modified LJ potential 



with energy and range parameters e m f and a w f, and a 5 w f for tuning the wetting property 
of the fluid. Both Uff and U w f are truncated at 2.5a. In our simulations, the density of 
fluid p equals to 0.81a~ 3 , the density of wall p w equals to 1.86a -3 (which determines the 
wall lattice constant), the mass of wall molecule m w equals to the mass of fluid molecule m, 
the parameters in U w f are e w f = 1.16e, a w f = 1.04a, and 5 w f = 1.0, and the temperature T 
is fixed at 2.8e/kB- The values of V w , L and H are varied as external conditions in different 
simulations. The steady-state flow fields are obtained from time averages over 10 4 to 10 6 r 



where r is the atomic time scale yma 2 /e. We have also performed similar simulations for 
other temperatures ranging from 1.2e//cs to 3.0e//ce- The MD velocity profiles can always 
be reproduced by our continuum model, with material parameters directly determined from 
MD simulations. This is due to the fact that the fluid remains to be Newtonian and the slip 
length is a constant at a given temperature. Details will be presented in Sec. VI B. 

We denote the region within z = 0.425a of the fluid surface the boundary layer (BL). 
It must be thin enough to ensure sufficient precision for measuring the slip velocity at the 
solid surface, but also thick enough to fully account for the tangential wall-fluid interaction 
force. The wall force can be singled out by separating the force on each fluid molecule into 
wall-fluid and fluid-fluid components. The fluid molecules in the BL, being close to the 
solid wall, can detect the discrete structure of the wall (the 'roughness' of the wall potential 
[3]). When coupled with kinetic collisions with the wall molecules, there arises a nonzero 






6 



tangential wall force density g™ that is sharply peaked at z ~ z$/2 and vanishes beyond 
Z ~ Zq . Here the subscript x in g™ and the z coordinates are for the BL at the lower 
fluid-solid interface (same below), with the understanding that the same physics holds at 
the other three fluid-solid interfaces. From the force density g™, we define the tangential 
wall force per unit area as G™(x) = J Z0 dzg™(x, z), which is the total tangential wall force 
accumulated across the BL. 

Spatial resolution along the x and z directions is achieved by evenly dividing the sampling 
region into bins, each Ax = 0.85(7 by Az = 0.425<r in size. The slip velocity v x ip in the 
lower/upper BL is obtained as the time average of fluid molecules' velocities in each BL bin, 
measured with respect to the moving wall (v* hp — v x +V w in the lower BL, or v* hp — v x — V w in 
the upper BL); the tangential wall force G™ in the lower/upper BL is obtained from the time 
average of the total tangential wall force experienced by the fluid molecules in each BL bin, 
divided by the bin area in the xy plane; the fluid stress component a xx ( zx ) is obtained from 
the time averages of the kinetic momentum transfer plus the fluid-fluid interaction forces 
across the constant-x(,2) bin surfaces, and the fluid velocity component v x ^ is measured 
as the time-average of that component within each bin. In particular, we have directly 
measured the fluid-fluid interaction forces across bin surfaces to obtain the contribution of 
intermolecular forces to the fluid stress, because the validity of the Irving-Kirkwood stress 
expression was noted to be not justified at a fluid-fluid or fluid-solid interface [12]. The 
technical details for the stress measurement near a fluid-fluid or fluid-solid interface may be 
found in Appendix B of Ref . [9] . 

III. SLIP BOUNDARY CONDITION 

Figure 2 shows the MD evidence for the existence of slip. It is seen that the slippage 
along x becomes quite large near the corner, as already observed in Ref. [10]. In particular, 
the fluid undergoes near-complete slip in approaching the corner, regardless of the system 
size. (By near-complete, we mean \vf ip \ approaches V w .) Far away from the corner (for fluid 
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that extends long enough along x, i.e., L 3> H), the flow is not perturbed by the vertical 
walls due to viscous damping, i.e., the uniform shear flow prevails and a small constant slip 
is detected. 

The NBC is the simplest alternative of the no-slip boundary condition. It states that 
the amount of slip is proportional to the tangential fluid stress at the solid surface. For a 
Newtonian fluid, the tangential viscous stress is proportional to the shear rate. Consequently, 
the NBC becomes that the amount of slip is proportional to the shear rate 7, i.e., v f ip = / s 7, 
where the proportionality constant l s is the slip length [3-5]. Physically, a nonzero slip length 
arises from the unequal wall and fluid densities, the weak wall-fluid interaction, and the high 
temperature. Together, they prevent the epitaxial locking of fluid layer(s) to the solid wall, 
and thus allow slip to occur. A recent study shows that fluid flow in carbon nanopores is 
characterized by a large slip length [13]. 

The verification of the NBC in the driven cavity flow, where slip velocity varies along 
the solid wall, consists of three stages, (i) We show that v s j %v is proportional to G™, the 
local tangential wall force per unit area, (ii) We show that G™ is balanced by the tangential 
fluid force G{.. This force balance is necessary because inertial effects are negligible in the 
BL of molecular thickness. Accordingly, v s J rp is also proportional to G{.. Here we emphasize 
that our stress measurement scheme has been designed to obtain the tangential fluid force 
correctly, (iii) We show that the shear stress is Newtonian. From (ii) and (iii) the slip length 
l s can be defined and measured. 

IV. MD RESULTS FOR HYDRODYNAMIC MODELING 

As shown in Fig. 3, the tangential wall force per unit area G™ is proportional to the 
local slip velocity v^ hp : 

g:{x) = -/*£*(*), (i) 

where the proportionality constant (3 is the slip coefficient. To see if nonlinearity would arise 
for large V w , MD simulations have been performed at V w ~ 1.0yje/m. Figure 3 shows that for 
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the wall-fluid interaction used here, nonlinearity in the slip boundary condition is accessible 
when \v x ip \ > 0.5^Je/m. According to the nonlinear effect discovered by Thompson and 
Troian [3], the slip length l s — r)/(3 increases with the increasing slip velocity when the latter 
is sufficiently large. This is seen in Fig. 3 where (3 decreases with the increasing \v s x hp \. 

The BL is thin enough to make the inertial term negligible in momentum equation 
(mpV w z /r] <C 1). It follows that the tangential wall force G x is balanced by the tangential 
fluid force per unit wall area, G[, i.e., G x (x) + G{.(x) = 0, as shown in Fig. 4. Here G{ is 
of the form 

f f z ° 
G s x (x) = a zx (x, zq) + d x / dza xx (x, z), (2) 

Jo 

coming from 

G l( x ) = I dz[d x a xx (x,z) +d z a zx (x,z)}, 
Jo 

together with the fact that a zx (x, 0) = 0. (More strictly, a zx (x,0~) = because there is no 
fluid below z — 0, hence no momentum transport across z = 0.) It follows from Eq. (1) and 
the tangential force balance that 

Gl(x) = 0v?*(x). (3) 

It is worth emphasizing that the normal stress a xx in the BL exhibits extremely large 
variation in the partial-slip region close to the corner, where v x hp varies quickly along x. In 
fact, the normal stress variation is so large that a small fluid density difference is even noted 
between the low- and high-pressure regions. Quantitatively, the BL-integrated normal stress 
Jq° dza xx (x, z) is essential to reaching Eq. (3), because close to the corner, the contribution 
of d x Jq° dza xx (x, z) to Gl is of the same order as that of a zx (x, z ). (In the region of uniform 
shear flow far away from the corner, d x Jq° dza xx (z) = 0.) 

To summarize, in obtaining the Navier slip condition Eq. (3), we first identify the BL and 
measure the tangential wall force therein to obtain Eq. (1). We then measure the normal 
and tangential stresses a xx and a zx according to the original definition of stress (because the 
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Irving-Kirkwood expression is not reliable near the fluid-solid interface). We finally calculate 
the tangential fluid force according to Eq. (2) to verify the equation of BL force balance. 

Equation (2) is due to the finite thickness of the BL, in which a tangential wall force is 
sharply distributed along the solid surface normal. Nevertheless, it is a reasonable expecta- 
tion that the fluid would experience almost the identical physical effect (s) from a wall force 
density G™6(z), concentrated strictly at z = with the same total wall force per unit area. 
Replacing a diffuse BL by a sharp BL can considerably simplify the form of the boundary 
condition, because local force balance along x then requires d x a xx + d z a zx = away from 
z = 0. Integration of this relation from + to z yields 



A comparison with Eq. (2) then relates G{ to a zx at surface: Gl(x) = cr zx (x, + ). Therefore, 
o zx changes from a zx (x, CT) = to a zx (x, + ) = G[(x) at z — 0, leading to 



This tangential fluid force density is in balance with the tangential wall force density G™6(z). 
Now the BL is from CT to + , instead of from to z as in the diffuse case. Correspondingly, 
the NBC becomes 



in the sharp boundary limit. 

It remains to be seen if the fluid is Newtonian, that is, if the viscous stress tensor is 
proportional to the rate of strain tensor. In particular, we need to find out if the tangential 
stress a zx is still proportional to the local shear rate 7 as the fluid-solid interface is ap- 
proached. We have measured both a zx and 7 at z = z , 2z , ■ ■ ■ (z = z Q is the top surface of 
the lower BL). As shown in Fig. 5, the ratio of a zx to 7 is indeed a constant at each z level, 
but this constant varies a little along z near the wall, due to the short-range density modu- 
lation induced by the rigid wall [14] . Far away from the wall, the fluid density approaches a 




(V-cr) -x 



G f J(z). 



a zx (x,0) = Pv s x hp (x) 



(4) 
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^-independent constant, and so does the ratio of a zx to 7. To see if non-Newtonian response 
would arise for large V w , MD simulations have been performed at V w ~ l.Oy/e/m. Fi Kure 
5 shows that the viscosity remains to be a constant for V w as large as \.2hyje/m. This is 
due to the fact that for V w ~ 1.0^/e/m and H > 10a, the shear rate is ~ O.lr -1 or smaller, 
whereas non-Newtonian response of bulk fluid is expected for 7 > 2t~ x . 

From Eq. (3) and the sharp boundary limit Gl(x) = a zx (x,0), we have Eq. (4). As 
the shear stress is verified to be Newtonian, i.e., a zx = 777, we obtain the commonly used 
hydrodynamic NBC 

vf p (x) = l s d z v x (x,0), (5) 

with the slip length l s = rj/(3. According to the sharp boundary limit involved in obtaining 
Eq. (5), we should compare the MD tangential fluid force G{.{x) with the continuum tan- 
gential viscous stress rjd z v x (x,0) at surface. Practically, a comparison between the MD and 
continuum profiles oiv* hp would suffice because G{ = (3v s J tp in the MD and r]d z v x (0) = (3v^ hp 
in the continuum hydrodynamics. 



V. CONTINUUM HYDRODYNAMIC MODEL 

A continuum hydrodynamic model has been formulated for the driven cavity flow from 
our knowledge of the NBC and the Newtonian stress. An explicit scheme has been designed 
to solve the hydrodynamic model, comprising the Navier-Stokes equation and the NBC. 
Material parameters include the fluid density p (= 0.81(T~ 3 ), the viscosity rj (= l.65y/em/a 2 ), 
and the slip length l s = r)/(3 (= 2.3a), all directly determined from MD simulations. Thus 
our model has no adjustable parameter. Numerical calculations show that steady-state flow 
fields from MD simulations can be quantitatively reproduced. 

The 2D flow is governed by the Navier-Stokes equation 



mp 



<9v 

- + (vV)v 



-Vp + 7]V 2 V, 
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with the incompressibility conditions V ■ v = 0, and the boundary conditions: v z — and 
lj 1 v* hp = —d n v x at the moving horizontal walls; v x — and l~ 1 vf w = —d n v z at the vertical 
walls (n denotes the outward surface normal). 

Our continuum model has six parameters, including the system dimensions L along x and 
H along z, the speed of the moving horizontal walls V w , the fluid density p, the viscosity rj, 
and the slip length l s . Taking H as the length unit, V w as the velocity unit, and i]V w /H as the 
pressure/stress unit, we are left with three dimensionless controlling parameters: the aspect 
ratio L/H, the dimensionless slip length l s /H, and the Reynolds number 1Z = pV w H/rj. In 
the regime of small Reynolds number (Stokes flow) the only controlling parameters are L/H 
and l s /H. 

The finite-difference scheme used for solving the Navier-Stokes equation is a modified 
version of the Pressure- Poisson formulation given in Ref. [15,9], where the incompressibil- 
ity condition is replaced by the pressure Poisson equation and a divergence-free boundary 
condition for the velocity. Variable grids with better resolution near the corners are used to 
save computational cost. 

VI. COMPARISON OF MD AND CONTINUUM RESULTS 

A. Stokes flow 

For 1Z ~ 1 or smaller, nonlinear effects associated with the inertial term are negligible. 
As a result, the dimensionless steady-state solution for the flow field, v(t/H)/V w , depends 
on L/H and l s /H only. Figure 2 shows the MD profiles of v x /V w in the BL (the slip profiles), 
obtained from three simulations using the same system size (L and H) but different wall 
speed V w . It is seen that the three MD profiles approach the same limiting profile [16]. This 
is due to the fact that the Reynolds number 1Z ranges from 0.067 to 8.3 in the three cases, 
thus justifies the Stokes-flow limit. In Fig. 6 the MD profiles of v x /V w at different z levels, 
obtained from two of the three simulations shown in Fig. 2, also indicate the Stokes-flow 
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limit. 

For comparison with the above MD results, the corresponding continuum results are also 
plotted in Figs. 2 and 6. They were calculated using the same set of material parameters 
p, T), and l s (see Sec. V) under respective conditions for L, H, and V w . The continuum 
calculations involve no adjustable parameter, and the overall agreement is satisfactory. A 
small discrepancy is noticed for the BL tangential velocity (or slip velocity) in a small region 
close to the corner where the slip amount is relatively large and displays sharp decay. This is 
presumably due to the short-range density modulation induced by the rigid wall [14], given 
the short distance H = 13.6a here [17]. In fact, this discrepancy tends to be less noticeable 
for larger H. 

B. Temperature effects 

MD simulations have been carried out as well for temperatures other than T = 2.8e//c#. 
We find that for T ranging from 1.2e/k B to 3.0e/k B , the MD velocity profiles can always 
be reproduced by our continuum model, with material parameters directly determined from 
MD simulations. In Fig. 7, we show the MD profiles of v x /V w at different z levels, ob- 
tained from a simulation at T = lAe/ks- The corresponding continuum results are also 
shown for comparison. They were calculated using the material parameters p = 0.81a -3 , 
T) = 1.7^/em/ a 2 , and l s = 1.27a. (MD measurements show that the viscosity weakly de- 
pends on the temperature, whereas the slip length is strongly temperature-dependent.) It 
is noticed that in the partial-slip region close to the corner, the discrepancy here is a bit 
larger than that seen for T = 2.8e//cs in Fig. 6. Again, this is caused by the short-range 
density modulation induced by the rigid wall [14], given the same short distance H = 13.6a 
here. In particular, such near-surface density modulation becomes more prominent as the 
temperature is lowered, and that's why the agreement here is less satisfactory. 

The uniform shear flow in the central region of the cavity is seen in Figs. 6 and 7. 
This part of the cavity flow is simply described by the Navier-Stokes equation d 2 v x = 
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plus the NBC v s ] %p = ±l s d z v x at z = and H, from which a constant slip amount Vq tp = 
2V w l s / (H + 2l s ) can be derived. Note this v S Q tp tends to vanish as H is sufficiently large. 

C. Power-law slip profile 

Now we turn to the variation of the slip velocity along the solid surface, still in the regime 
of small Reynolds number (the largest 1Z 13 for the largest H). We have performed a 
series of MD simulations using the same wall speed V w = 0.25^Je/m but different system 
size (L and H). The tangential slip velocity profiles at the fluid-solid interface, i.e., the slip 
profiles, are shown in the inset to Fig. 8. Regardless of the distance H, there is always a 
small core region in the immediate vicinity of the corner, on the order of a few l s , where the 
slip amount displays sharp decay. In particular, the structure of this core region is nearly 
independent of the system size. As H increases, however, a much slower variation of the slip 
profiles becomes apparent away from the core region. To find out the nature of this slow 
variation, we plot in Fig. 8 the same data in the log-log scale. The dashed line has the slope 
of —1, indicating a power-law behavior: away from the core region there is a wide partial-slip 
region in which the amount of slip varies as l/(x — x c ) where x c is the x coordinate of the 
corner. According to the NBC, the power-law variation of slippage means the same variation 
of tangential stress. Therefore, the asymptotic 1/r behavior of the stress variation [1,18] has 
indeed been observed in MD simulations. With H being finite and L being sufficiently large 
(L > 4if), far from the corners there is always a region of uniform shear flow, where the slip 
amount is a constant, given by v'q Hp = 2V w l s /(H + 21 s ). This is seen from the inset to Fig. 
8 where each slip profile shows a plateau. The small constant Vq %p (for H >> l s ) acts as an 
outer cutoff on the l/(x — x c ) profile at x — x c ~ H. For our largest MD simulation with 
H = 108.8cr (1Z m 13), the v shp oc l/(x — x c ) behavior actually extends to m 80a (or m 351 s ). 
Therefore, as H — > oo and Vq W approaches (no-slip), the power- law region can extend to 
hundreds of l s or even more. A large power-law partial-slip region is significant, because 
the outer cutoff length scale directly determines the integrated effects, such as the total 
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steady-state dissipation. While in the past the asymptotic 1/r stress variation away from 
the corner or the moving contact line has been known in continuum hydrodynamics [1,18], 
to our knowledge the observation that the partial slip is of the same spatial dependence has 
not been previously reported. 

The continuum results are also shown in Fig. 8 for comparison. They were calculated 
using the same set of material parameters p, rj, and l s corresponding to the same local 
properties in all the four MD simulations. The overall agreement with the MD results is 
excellent. This not only clearly demonstrates the validity of our continuum model, but also 
confirms the power-law partial-slip region in the continuum hydrodynamics. 



MD simulations have been carried out to investigate the flow fields at large Reynolds 
number. Figure 9 shows the MD profiles of v x /V w at different z levels, obtained from a large- 



scale simulation for V w = 0.5ye/m, L = 511(7, and H = 204cr (TZ 50). The approximate 
fore-aft symmetry of the Stokes flow disappears. In particular, the slip profiles associated 
with the left and right corners are no longer symmetric, especially in the partial-slip regions 
where the 1/r variation would appear were the Reynolds number being small. This deviation 
from the power-law behavior described in Sec. VI C is due to the vorticity convected with 
the fluid. The continuum profiles of v x /V w are also shown in Fig. 9. Excellent agreement 
is seen, from the fine features in the molecular-scale vicinity of the solid walls to the global 
flow. 



In summary, we have carried out MD simulations to study the fine structure of the driven 
cavity flow. It has been verified that the Navier boundary condition can quantitatively 
describe the fluid slipping at the solid wall. It has also been shown that close to the corner, 
where velocity variation is extremely fast, the shear stress is still Newtonian. Based on these 



D. Large Reynolds number 




VII. CONCLUDING REMARKS 
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MD facts, a continuum hydrodynamic model has been formulated. This model involves no 
adjustable parameter, and it can produce flow fields in quantitative agreement with those 
from MD simulations. 

Recently, people have developed some MD-continuum hybrid methods for the study 
of fluid dynamics that involves complex small-scale structure where validity of continuum 
formulations is not clear [19-22]. These hybrid methods have been successfully applied 
to study the moving contact line [20], channel flow with nano-scale rough wall [21], and 
the corner singularity in driven cavity flow [22]. Here we want to point out that when 
formulated correctly, continuum hydrodynamic approach may still be applicable to some 
of these problems, including the moving contact line [9] and the driven cavity flow. For 
each problem, the validity of its continuum model has been verified, first by a direct MD 
measurement and then by a comparison with full MD results. 
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FIG. 1. Geometry of MD simulations. The empty circles indicate the instantaneous molecular 
positions of the fluid projected onto the xz plane. The solid squares denote the wall molecules. 
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FIG. 2. Boundary-layer tangential velocity profiles. The scaled tangential velocity v x /V w in the 
lower BL is plotted as a function of x/a. (That for the upper BL can be obtained by symmetry.) 
The lower wall is moving at —V w , hence v x /V w = means complete slip and v x /V w = —1 means 
no slip. Close to the corner (within a few a) the slip is near-complete while away from the corner 
the slip becomes smaller. There is a uniform shear flow in the central region where the slip is a 
constant. The three cases shown here are of the same L = 61.3a and H = 13.6cr but different V w . 
The circles, squares, and diamonds denote the MD results for V w = 0.25i/e/m, 0.01\/e/m, and 
1.25-y/e/m, respectively. (The squares for the smallest V w show the largest statistical fluctuation.) 
It is seen that the three MD profiles show negligible deviation from each other. The solid, dashed, 
and dotted lines denote the continuum results for V w = 0.25ye/m, O.Olye/m, and 1.25-y/ e/m, 
respectively. They also show negligible deviation from each other. 
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FIG. 3. The tangential wall force G™ (in the unit of e/cr 3 ) is plotted as a function of the slip 
velocity v^ ip (in the unit of y/e/m) for the lower {y^ ip > 0) and upper (vf p < 0) BL's. The two 
cases shown here are of the same L = 61.3a and H = 13.6a but different V w . The circles and 



diamonds denote the cases of V w = 0.25-^/e/m and 1.2by/e/m, respectively. The dashed line is for 
eye guidance. Linearity is seen to be well preserved for V w = 0.25^ e/m (circles, see also the inset 
for an enlarged plot). The slip coefficient (3 is obtained to be 0.69^/em/a 3 . For V w = 1.25\/e/m 
(diamonds), however, nonlinearity shows up for \v^ tp \ > 0.5y/e/m. 
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FIG. 4. The tangential wall force G™ (in the unit of e/cr 3 ) is plotted as a function of the 
tangential fluid force G{. (in the unit of e/cr 3 ) for the lower (G[ > 0) and upper (G£ < 0) BL's. 
The two cases shown here are of the same L = 61.3a and H = 13.6a but different V w . The circles 



and diamonds denote the cases of V w = 0.25^/e/m and 1.25-\/e/m, respectively. The dashed line 
has the slope of —1, indicating G™ + Gl = 0. The inset shows an enlarged plot for V w = Q.2b^J7Jm. 
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FIG. 5. Newtonian response of shear viscous stress. The tangential stress a zx (in the unit of 
e/cr 3 ) is plotted as a function of the shear rate d z v x + d x v z (in the unit of t -1 ). The upper panel 



shows the case of L = 61.3a, H = 13.6(7, and V w = Q.2b^Je/m] the lower panel shows the case of 



L = 61.3<7, H = 13.6o", and V w = 1.25^/ e/m. In each panel, the circles denote the data collected at 
the levels of z = zq and z = H — zq, the squares denote the data collected at the levels of z = 2zq 
and z = H — 2zo, the diamonds denote the data collected at the levels of z = 3zo and z = H — 3zo, 
the triangles denote the data collected at the levels of z = Azq and z = H — Azq, and the pluses 
denote the data collected at all other levels. At each z level, a zx shows a linear dependence on 
d z v x + d x v z . The ratio of a zx to d z v x + d x v z varies across the first four z levels away from the 
wall, and becomes a z-independent constant at levels deeper in the fluid. Note that the ratio of 



a zx to d z v x + d x v z remains to be unchanged from V w = 0.25^/ e/m to V w = 1.25^/ e/m. (The scales 
for a zx and d z v x + d x v z in the lower panel are five times larger than those in the upper panel, 
corresponding to the five times difference between the two values of V w .) 
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FIG. 6. Comparison between the MD and continuum results for two Stokes flows. The sym- 
bols denote the MD profiles of v x /V w at different z levels, obtained for the same L = 61.3a and 
H = 13.6<t but different V w . The empty and solid symbols represent the cases of V w = 0.25\/e/m 
and 0.01y/e/m, respectively. The v x /V w profiles are symmetric about the center plane z = H/2, 
hence only the lower half is shown for z = 0.2125<7 (circles), 1. 9125c (squares), 3.6125<r (diamonds), 
and 5.3125(7 (triangles). Large fluctuation is seen from the solid symbols for V w = O.Olye/m 
(smaller than the thermal velocity by two orders of magnitude), although we averaged over ~ 10 6 r 
to reduce noise. The solid and dashed lines denote the corresponding continuum profiles for 
V w = 0.25VeM = 1-67) and 0.01 ^/IJm (K = 0.067), respectively. 
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FIG. 7. Comparison between the MD and continuum results for a Stokes flow at temperature 
T = lAe/ks- The symbols denote the MD profiles of v x /V w at different z levels, obtained for 
V w = 0.25ye/m, L = 61.3a, and H = 13.6<r. The v x /V w profiles are symmetric about the center 
plane z = H/2, hence only the lower half is shown for z = 0.2125<t (circles), 1.9125cr (squares), 
3.6125cr (diamonds), and 5.3125cr (triangles). The solid lines denote the corresponding continuum 
profiles. 
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FIG. 8. Log-log plot of the slip profiles showing the power-law behavior. Here v x /V w + 1 is the 
scaled slip velocity in the lower BL, and {x — x c )/a measures the distance from the corner, with 
x c = —L/2 being the coordinate of the left corner and x ranging from x c to (the cavity center). 
The lower wall is moving at —V w , hence v x /V w +l = 1 means complete slip and v x /V w +l — > means 
no slip. The MD results were obtained from four simulations for the same V w = 0.25-\/e/m, but 
different H. The symbols represent the MD results and the lines represent the continuum results, 
obtained for H = 13.6<t (circles and solid line), H = 27.2a (squares and dashed line), H = 54. 4a 
(diamonds and dotted line), H = 108. 8cr (triangles and dash-dotted line). The thin straight line 
has the slope of —1, indicating that the l/(x — x c ) behavior is approached for increasingly larger 
H. For H = IO8.80", the power-law behavior extends from x — x c ~ 14a ~ 6l s to 80<r ~ 35l s . Inset: 
The scaled tangential velocity v x /V w in the lower BL, plotted as a function of (x — x c )/a. 
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FIG. 9. Comparison between the MD and continuum results for a flow at TZ ~ 50. The symbols 
denote the MD profiles of v x /V w at different z levels, obtained for V w = 0.5\/e/m, L = 511<r, and 
H = 204<7. The v x /V w profiles are symmetric about the center plane z = H/2, hence only 
the lower half is shown for z = 0.425<7 (circles), 17.425<r (squares), 34.425<7 (diamonds), 51.425<r 
(up triangles), 68.425<r (down triangles), and 85.425a (left triangles). The solid lines denote the 
corresponding continuum profiles. 
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